% This file generates Figure B1 in Appendix B of the paper. It presents the
% IRF as a function of the duration of ZLB from a stochastic peg.

close all; clc; clear;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Parameters and Steady-State Values %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

beta                    =   0.99;   % Discount factor.
chi                     =   0;      % Inverse of Frisch elasticity. 
alpha                   =   0.8;    % Probability interest rate is fixed.
sigma                   =   1;      % Inverse of EIS.

phi                     =   0.75;   % Prob. firms cannot adjust prices. 

rho                     =   0.9;    % Persistence of natural rate shock. 

% Slope Phillips Curve:
gamma                   =  (1-phi)*(1-phi*beta)*(sigma + chi)/phi;

% Output gap policy:
theta2                  =   (1-alpha*beta*rho)/(sigma*(1-alpha*rho)*(1-alpha*beta*rho)...
                             - alpha*rho*gamma);
% Output gap policy:
theta1                  =   gamma/(sigma*(1-alpha*rho)*(1-alpha*beta*rho)...
                            - alpha*rho*gamma); 

H                       =   21+1;     % Number of periods.

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Creating IRFs                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Pre-allocating:
rf_irf                  =   zeros(1,H);
x_irf                   =   zeros(1,H);
y_irf                   =   zeros(1,H);
A_irf                   =   zeros(1,H);
yf_irf                  =   zeros(1,H);
pi_irf                  =   zeros(1,H);

% Initial values:
rf_irf(1)               =   sigma*(1+chi)*(rho-1)/(sigma +chi);
pi_irf(1)               =   theta1*rf_irf(1);
x_irf(1)                =   theta2*rf_irf(1);
A_irf(1)                =   1;
y_irf(1)                =   x_irf(1) + (1+chi)/(sigma+chi)*A_irf(1);
yf_irf(1)               =   (1+chi)/(sigma+chi)*A_irf(1);

% Creating IRFs:
for jx = 2:H
    rf_irf(jx)              =   rho*rf_irf(jx-1);
    A_irf(jx)               =   rho*A_irf(jx-1);
    pi_irf(jx)              =   alpha^(jx-1)*theta1*rf_irf(jx);
    x_irf(jx)               =   alpha^(jx-1)*theta2*rf_irf(jx);
    y_irf(jx)               =   x_irf(jx) + (1+chi)/(sigma+chi)*A_irf(jx);
    yf_irf(jx)              =   (1+chi)/(sigma+chi)*A_irf(jx);
end

%%
% Considering a different value of Alpha. 
% Now: 
alpha                   =   2/3;
theta2                  =   (1-alpha*beta*rho)/(sigma*(1-alpha*rho)*(1-alpha*beta*rho)...
                            -alpha*rho*gamma); 
theta1                  =   gamma/(sigma*(1-alpha*rho)*(1-alpha*beta*rho)-alpha*rho*gamma); 

% Pre-allocating:

pi_irf2                 =   zeros(1,H);
x_irf2                  =   zeros(1,H);
y_irf2                  =   zeros(1,H);
A_irf2                  =   zeros(1,H);
yf_irf2                 =   zeros(1,H);
% Initial values
pi_irf2(1)              =   theta1*rf_irf(1);
x_irf2(1)               =   theta2*rf_irf(1);
A_irf2(1)               =   1;
y_irf2(1)               =   x_irf2(1) + (1+chi)/(sigma+chi)*A_irf2(1);
yf_irf2(1)              =   (1+chi)/(sigma+chi)*A_irf2(1);

% Creating IRFs:
for jx = 2:H
    A_irf2(jx)              =   rho*A_irf2(jx-1);
    pi_irf2(jx)             =   alpha^(jx-1)*theta1*rf_irf(jx);
    x_irf2(jx)              =   alpha^(jx-1)*theta2*rf_irf(jx);
    y_irf2(jx)              =   x_irf2(jx) + (1+chi)/(sigma+chi)*A_irf2(jx);
    yf_irf2(jx)             =   (1+chi)/(sigma+chi)*A_irf2(jx);
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Plotting results                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t                       =   0:H-2;

figure
subplot(1,2,1)
plot(t,yf_irf(1:H-1),'-b',t,y_irf2(1:H-1),'--b',t,y_irf(1:H-1),':b','Linewidth',1.5)
xlabel('Horizon','Interpreter','latex','fontSize',14)
ylabel('Output','Interpreter','latex','fontSize',14)
h = legend('No ZLB','p = 2/3 ','p = 4/5','Location','Best');
    set(h,'Interpreter','latex')
grid on
set(gca,'layer','top')

subplot(1,2,2)
plot(t,zeros(1,H-1),'-b',t,pi_irf2(2:H),'--b',t,pi_irf(2:H),':b','Linewidth',1.5)
axis([0 H-2 min(pi_irf)-0.1 0.1 ])
xlabel('Horizon','Interpreter','latex','fontSize',14)
ylabel('Expected Inflation','Interpreter','latex','fontSize',14)
h = legend('No ZLB','p = 2/3 ','p = 4/5','Location','Best');
    set(h,'Interpreter','latex')
grid on
set(gca,'layer','top')


